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This work provides a ground for a quantitative interpretation of experiments on step bunching 
during sublimation of crystals with a pronounced Ehrlich-Schwoebel (ES) barrier in the regime 
of weak desorption. A strong step bunching instability takes place when the kinetic length dd = 
D s /Kd is larger than the average distance I between the steps on the vicinal surface; here D s is 
the surface diffusion coefficient and Kd is the step kinetic coefficient. In the opposite limit dd <C I 
the instability is weak and step bunching can occur only when the magnitude of step-step repulsion 
is small. The central result are power law relations of the form L ~ H a , l m i n ~ H~~" between 
the width L, the height H, and the minimum interstep distance l m in of a bunch. These relations 
are obtained from a continuum evolution equation for the surface profile, which is derived from the 
discrete step dynamical equations for the case dd 3> I- The analysis of the continuum equation 
reveals the existence of two types of stationary bunch profiles with different scaling properties. 
Through comparison with numerical simulations of the discrete step equations, we establish the 
value 7 = 2/(n + 1) for the scaling exponent of i m i n in terms of the exponent n of the repulsive 
step-step interaction, and provide an exact expression for the prefactor in terms of the energetic and 
kinetic parameters of the system. For the bunch width L we observe significant deviations from the 
expected scaling with exponent 7 = 1 — 1/a, which are attributed to the pronounced asymmetry 
between the leading and the trailing edges of the bunch, and the fact that bunches move. Through 
a mathematical equivalence on the level of the discrete step equations as well as on the continuum 
level, our results carry over to the problems of step bunching induced by growth with a strong inverse 
ES effect, and by electromigration in the attachment/detachment limited regime. Thus our work 
provides support for the existence of universality classes of step bunching instabilities [A. Pimpinelli 
et al., Phys. Rev. Lett. 88, 206103 (2002)], but some aspects of the universality scenario need to 
be revised. 

PACS numbers: 68.35.-p, 81.10.-h, 05.70.Np, 89.75.Da 

I. INTRODUCTION 

The formation of step bunches at a vicinal surface is a problem of great current interest, both from a fundamental 
viewpoint and with regard to the possible uses of step bunches as nanotemplates or nanostructuresi2i2iiSi&. Mech- 
anisms causing step bunching instabilities include strain effectsi*2*2i£, sublimation under conditions of asymmetric 
detachment kinetics known as the Ehrlich-Schwoebel (ES) effecfe^iiSiii, growth with an inverse ES effec1i2ii24iii4 I and 
surface electromigration 15 i 16 . 17 i 18 i 19 i 20 . 21 i 22 i 23 i 24 i 25 i 26 i 27 i 28 . 

Quite recently it was realised that step bunching is a promising way to study the interactions between the 
steps 29 ' 30 ' 31 ' 32 ' 33 ' 34 . The physical ground is simple: The steps in the bunch keep a certain distance from each other 
because the step-step repulsion balances the tendency to further compression of the bunch. The free energy related 
to the step-step interaction is of the form A/l n , where I is the interstep distance. When n — 2, the amplitude A(T) 
accounts for both elastic and entropic repulsion between the steps 3 *. Under crystal-vapour equilibrium one has the 
relation A(T) ~ g(T) where g(T) is the step repulsion coefficient in the expression 

/(p) = f(Q) + K P + gp 3 (1) 

for the surface free energy (per unit projected area) of a vicinal crystal surface with a density of steps p. 

To infer information about the step-step interactions from experimental observations of bunch morphology, one 
makes use of scaling relations between the length and time scales characterizing the bunches. The relevant length 
scales are the width L and the height H of the bunch, and the spacing £ between subsequent bunches (FigJIJ. The 
length £ is also sometimes referred to as the terrace width; this nomenclature is somewhat ambiguous, however, 
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FIG. 1: Schematic of a bunched vicinal surface, illustrating the definition of the bunch width L, the bunch height H and the 
bunch spacing £. 



because the region between two bunches may contain several monoatomic steps, and, hence, several wide terraces. 
The bunch height is related to the number of steps N in the bunch by H = Nho, where ho is the height of an atomic 
step. A quantity that is directly accessible to experimental observations^ is the minimal terrace size i m i n inside the 
bunch, which is related to the maximal slope m max through l m i n = /jo/m max . Following the notation of Ref.36, we 
introduce scaling exponents a and 7 characterizing the shape of individual bunches through the relations 

H~L a , U~iV-l (2) 

Assuming that the minimal terrace size Z m ; n is of the same order as the mean size I = L /N of the terraces inside the 
bunch, one expects the exponent identity 

7=1- l/a. (3) 
Furthermore the coarsening of the bunch morphology with time is described by a dynamic exponent z defined through^ 6 . 



£~N~t a / z . (4) 

Because the ratio £/N has to be equal to the mean terrace width I, which is fixed by the overall vicinality of the 
surface, the bunch spacing £ and the bunch size N grow with time in the same manner. 

For step bunching induced by surface electromigration, scaling relations of the form J5J have been derived both 
for non-transparent and transparent steps^^SiSL^. Their application to the analysis of experiments^i^ proceeds 
in two stages. First, the value of the scaling exponent a or 7 is used to determine the kinetic regime and the value 
of the step interaction exponent n, and subsequently an estimate for the strength of the step-step interaction (in 
relation to the driving force for step bunching) is extracted from the prefactor of the scaling relation. For example, for 
non-transparent steps one finds that a = (n + l)/(n— 1), which yields a = 3 in the standard case n = 2, whereas for 
n = 1 the width of the bunch L would in fact be independent of the bunch size. Provided that n = 2 and a = 3, the 
bunch width depends on the step interaction coefficient as L ~ LA(T)] 1 / 3 . This relation provides a ground to study 
the temperature dependence of A(T) by measuring the width of the bunch as a function of temperature^. 

The main purpose of this paper to derive scaling relations of the form for step bunches induced by the ES 
effect during sublimation, and to put the results into perspective with regard to previous work on other step bunching 
instabilities. In the next section we introduce the basic concepts and quantities of the Burton, Cabrera, Frank (BCF) 
models in the presence of ES barriers (see Ref.38 for a recent review). In Sect IIIII the equations of step motion are 



3 



displayed and various limiting cases are discussed. In Sections IIVI and |V| a continuum evolution equation for the 
surface is derived, and the structure of stationary bunch solutions is analysed. A key result is the existence of two 
types of solutions with different scaling properties in the sense of (0). In Sect IVII the mathematical equivalence of 
the sublimation problem to appropriate limiting cases of step bunching induced by electromigration and growth with 
inverse ES barriers is pointed out and exploited. Sectior fVlII presents results obtained from numerical simulations of 
the discrete step dynamics and compares them to the predictions of the continuum theory. In Sect lVUll we critically 
examine 3 ^ a recently proposed classification scheme for step bunching instabilities 36 in the light of our results, and 
some general conclusions are drawn in Sect llXI 



II. DISCRETE MODEL AND BASIC CONCEPTS OF SUBLIMATION BY STEP FLOW 

We consider a vicinal surface going uphill in the +x direction. The processes of atom migration (in the presence of 
desorption and deposition) are described by the stationary diffusion equation 

„ d 2 rie n, „ 

Ds-r^- — + i? = (5) 

ax t s 

where D s is the coefficient of surface diffusion, n s is the concentration of mobile atoms, adsorbed on the surface, R is 
the deposition rate of atoms to the crystal surface, r s is the life time of an atom in a state of mobile adsorption, and 
x is a coordinate perpendicular to the step edges (we consider a system of parallel steps with straight edges). The 
exchange of atoms between the crystal phase and the dilute layer of adatoms takes place at the steps and determines 
the boundary conditions for Eq. (JSJ). It is essential to note that the ES barrier considerably decreases the permeability 
of the steps. Really, an atom has to break many chemical bonds when it crosses over the step down to a position of 
adsorption at the step edge. Therefore, such an event happens very rarely. The opposite jump from the position of 
adsorption at the step edge to a position of adsorption on the upper terrace (or, briefly, from the lower to the upper 
terrace) also happens rarely in accordance to the principle of the detailed balance. This circumstance justifies the 
reduction of the diffusion problem on the crystal surface to a diffusion problem at a single terrace. The boundary 
conditions relate the surface fluxes of adatoms with the power of the step as a generator (or a sink) of adatoms. For 
the terrace between the i-th and i + I-th step one can write 

dfi 

-Ds—j^-\ x =xi = ~K u [n s {xi) - n e s (xi)] (6) 

-D s ^-\ x = Xi+1 = K d [n s (x i+ i) - n e s (x i+ x)] 

where K ( i and K u are the kinetic coefficients for ascending and descending steps. The step kinetic coefficients and 
K u are defined by the expression for the velocity of the step motion 

dx ■ 

v= -rj- = -a±a\\{K d [n St i-i(xi) - n% (Xi)] + K u [n s<i (xi) - n e s (x l )}} (7) 

where a± and a|| are the interatomic distances perpendicular and parallel to the steps, n St i-i(x) is the adatom 
concentration on the terrace between the i — 1-th and the z-th step, and n e s (xi) is the equilibrium value of the adatom 
concentration in the vicinity of the step, situated at x\. The value of n e s in the vicinity of the i-th step depends on 
the distances to the neighboring steps i + 1 and i — 1, since n e s (xi) = n e s exp[Afi(xi)/kBT], wher e - 29 ; 30 



Afi(xi) ( l \ ( l 



kj$T \Xi+i Xi J \Xi Xi—i 

and 



(8) 



/204(T)\ 



1/3 
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is a length scale characterizing the strength of the step-step interaction. In © we have introduced the atomic area 

fi = <2j_<Z|| . 

The Ehrlich-Schwoebel harried provides the physical ground for the inequality K u ^ Kd, since the exchange of 
atoms between the crystal phase and the adlayer on the upper terrace has a lower rate than the exchange between the 
crystal and the adlayer at the lower terrace, i.e., K u < Kd, although the opposite inequality has also been discussed^. 
Assuming that the step kinetic coefficients can be written in the form K v ,d — Koexp(—E Ui d/kBT) one can write 
K u /Kd = exp(— AE/k^T) where AE = E u — Ed characterises the asymmetry in the atom attachment-detachment 
kinetics at the steps. The ratio (3 = K u /Kd is an essential parameter in the considerations presented below. 



III. EQUATIONS OF STEP MOTION 



Since the normal ES barrier (AE — E u — Ed>0) causes a step bunching instability only in the case of sublimation 9 , 
the equations of step motion will be derived under the condition R = (i. e., in the absence of deposition). The 
steps then move to the right. Solving the diffusion problem [Eq. JSJ) with the boundary conditions JBJl] one obtains 
an expression for the adatom concentration at the crystal surface. Substituting n s j-i(xi) and n St i(xi) into Eq. (JZJ) 
one can write an expression for the velocity dxi/dt of the i-th step as a function of the widths of the lower (left) and 
upper (right) terraces 



dxi 
~dt 



Vd- 



(10) 



In the physically interesting limit of small desorption rate, in the sense that the diffusion length A s = \/D s t s 3> 
Xi+\ — Xi, the two contributions to the step velocity l|10(l read 



£(l + /3) + 



and 



I'd 



k-aT k B T 



(ii) 



(12) 



Here we have introduced the kinetic length; 
(|11I12|1 have two limiting cases: 



miLdd 



(a) ^(l+/3) « 



D s /Kd and d u = D s /K u = D s /Kd(3 = dd/(3. The equations 

(13) 



A, 



(b) £(1 

As 



/?)» 



Xi-l 



A., 



(14) 



The limit (a) (realised when Xi — £&j_i ^> dd and dd/X s <C 1) reduces the denominators of Eqs. (|11I12|> to ((3/X s )(xi — 
Xi-i). If in addition {3 3> (dd/A s ) 2 , the parameter (3 is eliminated from the terrace-length dependent destabilizing 
terms in (|llf) and (|12f) . and appears only in a constant contribution to the step velocity, where it does not provide 
any ground for a step bunching instability. It is, however, quite possible that an instability is induced by higher order 
terms [(xi — a;i_i)/A s ]' y with v > 2. We shall address this possibility in Scct lVIII devoted to numerical analysis of the 
discrete model. 

The limit (b) is more interesting. It takes place under the assumption x% — Xi-i <C dd and dd/X s < 1. Then only 
the constant terms in the denominator of Eas. (|llll2fl are retained and the terms quadratic in the terrace lengths can 
be neglected relative to the linear ones. Thus Eqs. H11I12|I reduce to 



_ D a Sln e a f P 



d d {l+f3) \k B T 



[Afi(xi) - A^(x l+ i) 



(3d d , 

-^2-(s<+l 



Xi) 



(15) 
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Vd 



d d (l+P) \k B T 



(16) 



The two terms in the curly brackets have a clear physical meaning. The first term (the difference between the chemical 
potentials of neighbouring steps) is the driving force for the relaxation of the fluctuations in the step density. The 
second term in Eqs. 1)15116(1 reflects the asymmetry in the step kinetics and provides a ground for step bunching 
instability. These two terms describe a model previously analysed in Ref.42, where it was shown how continuum 
equations for the one-dimensional crystal profile h(x, t) can be derived exactly, provided the step dynamics is linear 
in the terrace lengths. Indeed, apart from the chemical potential difference in the square brackets, we can write 



dxi 
~dt 



V u + V d 



D a Q,n% 



P , 1 / 



where 



D s Sln e s Vtn e s _ - 

= = ilR e = R e 

A? To 



(17) 



(18) 



is the equilibrium value of the desorption rate per adsorption site. We will return to the continuum equation derived 
from below in Sect El 

IV. CONTINUUM LIMIT OF THE RELAXATIONAL DYNAMICS 

In this section we develop a continuum description for the relaxational part of the step dynamics. The chemical 
potential differences in Eqs. i|15H6fl can be written approximately in the form 



1 



[Afj,(xi) - An(xi-i)] 



1 



d_ 

dx 



k B T 1 ^ v ' ^ " k B T 

Making use of the relation (xi — Xi-i) rs ho{dh / dx)^ 1 , one can bring the last expression into the form 

-l 



ho ( dh 



k B T \ dx 



d_ 

dx 



Afifa-i] 



(19) 



(20) 



This is the contribution to the rate of motion of the i-th step, due to the difference between the chemical potentials 
of the i-th and the i — 1-th step. In the same way the contribution of the difference between the chemical potentials 
of the i-th and i + 1-th step to the motion of the i-th. step is 



hp 

k B T \dx ) 



dhX 



(21) 



Thus the total contribution of the variation of the chemical potential to the rate of step motion is obtained by 
substituting Eqs. (EOJl and QU into Eqs. I|15I16|I and (fT77|) . 



dxi D s flho n e a (3(xi — d J f dh 

~dT ~ dd(l + 0)k B T dx I \ dx~ 



d_ 

dx 



Afi(xi 



(22) 



It can be seen from Eq. 1)221) that the rate of relaxation of a non-equilibrium configuration of steps towards an ideal 
equidistant configuration (having a zero contribution of the step-step repulsion to the chemical potential) is low when 
the parameter (3 is small, i.e., when the ES barrier is high. This is easy to understand since the relaxation of the step 
configuration takes place by detachment of atoms from those steps with high chemical potential, and their subsequent 
attachment to other steps with low chemical potential. These atoms have to overcome the ES barrier either in the 
detachment, or in the attachment process. 

Substituting the expression (12211 into the equation 



dh 
~dt 



h 



dxi I dt 



(23) 
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used by Frank - in developing the kinematic theory of crystal growth, one obtains 



dh 



d_ 

dx 



where 



d d (l + [3)k B T 



Oh 



(24) 



(25) 



is a surface mobility relating the mass current in the curly brackets on the right hand side of i|24|) to the gradient of 
the chemical potential A/x. The proportionality of the mobility (|25|) to the inverse of the surface slope is intimately 
related 4 ^ 5 - to our assumption of a large kinetic length (slow detachment/attachment), i.e., case (b) [Ea. (|14fl ] . Indeed, 
carrying out the same manipulations for the opposite case (a) [Eq. l|13|) ]. one arrives instead at the expression 



D s nh n e s (3 



k B T[[3 + (d d /X s 



(26) 



which is independent of the surface slope. 

Since the continuum limit of the expression JHJ) is^i 



Eq. lt^|) can be presented as 



A/z 



dh + d_ \ c (9h 
dt dx | \ dx 



6nA(T) dh d 2 h 
h 2 , dx dx 2 ' 



d_ 

dx 



dh d 2 h 
dx dx 2 



= 



(27) 



(28) 



where 



C 



6A(T)D s n 2 n% (3 _ 3l%nD s n e J 
d d {\+l3)k B T ~ d d (l + f3) ' 



(29) 



A similar equation for the relaxation of the step bunches has been published previously 45 . It describes the relaxation 
of the fluctuations in the step density due to the step-step repulsion. Equation l|28|l is highly nonlinear and differs 
qualitatively from the linear evolution equation originally introduced by Mullins^ to describe the relaxation of a 
crystal surface above the roughening temperature. 



V. CONTINUUM THEORY OF BUNCH SHAPES 



A. The continuum evolution equation 

It was shown in Ref.42 how the continuum limit for a set of step equations of the linear form (|17|l can be carried out 
in an essentially rigorous manner. The main idea is to perform the coarse graining operation on the level of the linear 
step dynamics, where it can be done exactly, and subsequently derive the evolution equation for the height profile 
h(x, t) through a nonlinear variable transformation. Combining the result of this procedure with the continuum limit 
(|28|l for the relaxational dynamics, we obtain the full evolution equation for our problem. It takes the form of a 
continuity equation, 

dh dJ . 

at + ITx - - hoR °> (30) 

where the expression for the current is 

_ R e h 2 (l - 0) 1 RJ4 1 dm CI d 2 m2 

2(1 + 0) m 6 m 3 dx 2 m dx 2 ' 1 ' 

Here m = dh/dx is the surface slope, which will be taken to be positive. The first term on the right hand side of 
(|31|l is the destabilizing, downhill current which is responsible for the step bunching instability, while the last term 
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describes the smoothening effect of the step-step repulsion, as described above. The second term is the only one to 
break the reflection symmetry (x — ► —x and m — * —to) of the evolution equation. As will be shown later, this term 
leads to different behavior at the upper and the lower edges of the bunch. Following earlier work^ 2 *^, we will refer to 
it as the symmetry-breaking term. 

It is noteworthy that the continuum equation l|30|l conserves the volume of the him, apart from the constant 
sublimation rate on the right hand side, which does not couple to the surface morphology. A dependence of the 
sublimation rate on the surface slope will be felt when the step spacing becomes comparable to the diffusion length^l. 
In this sense Eq. l|5U|l is valid only to leading order in l/X s . 



B. Linear stability analysis 

It is straightforward to derive from (|30I31|1 the instability condition for a vicinal surface of slope mo = ho /I. We set 
h(x, t) = m a x — R e h a t + e q (x, t), where e q (x, t) ~ e iqx+u(q)t j g a p er t U rbation of wavenumber q, and expand (|3UI31|) 
to linear order in e. This yields the expression 



Re.h (l (3) 2 4 R e hl 3 

^ 2(1 + /3K 9 - Cq ~ 6^T 1<? (32) 
for the growth rate of the perturbation. The perturbation grows when 1Z(lu) > 0, i.e. for wavenumbers q < <7 max = 



Reh^l — /3) / [2(1 + /3)toqC], and perturbations with wavenumber q* = q max /V2 are maximally amplified. The 
corresponding wavelength is 

A* = ^ = 4nV3Sl, (33) 

where we have introduced the dimensionless quantity 

g Ag Zg K d K u r fl \ 3 
1 - (3 ckl* K d -K u l\l ) ■ 1 ' 

We will see below that the physical parameters of the problem enter the properties of the bunch shape only in 
the combination (|34|l . It is interesting to note that S does not depend on the surface diffusion constant D s . The 
wavelength A* determines the linear size of bunches in the beginning of the bunching instability. Correspondingly the 
number of steps in an incipient bunch is given by 

N* = X*/l = AttV3S w 21.8 x Vs. (35) 

The imaginary part of the growth rate l|32|) . which derives from the symmetry-breaking term in l|31|) . does not affect 
the stability of the surface, but it induces a drift of fluctuations. 



C. The mechanical analog 

Our main goal in this section is to compute the shape of large, almost stationary bunches from Ea. l|31|) . For a 
stationary profile, the current jSJ is set to a constant Jo. We neglect the symmetry-breaking term for now, and 
return to its relevance at the end of the section. Introducing the quantity u = m 2 , which is positive by construction, 
the stationarity condition J = Jo can then be written in the familiar form 

Cd 2 u 1/9 dV 

of a classical particle coordinate u{x) moving in a potential V(u) = -Bu- (2/3) J u 3 / 2 , where B = Refill- j3)/ {2(1 + 
(3)] > and Jo has to be chosen negative. The bunch is a particle trajectory which starts at u — at "time" x = 0, 
reaches a turning point u = u max at time x = L/2, and returns to u — at time x = L. Here L can be identified with 
the bunch width, and the minimal terrace length in the bunch introduced previously is given by ^ m i n = ho/y/u max . 

As can be seen in Fig. the potential V(u) admits two types of periodic trajectories. When the total particle 
"energy" £ = (C /4)(du/dx) 2 + V(u) is negative, u performs an essentially harmonic motion around the minimum of 
the potential. This will be referred to as a type I solution. The particle velocity du/dx vanishes at both turning points, 




FIG. 2: Sketch of the potential V(u) appearing in the mechanical analog. The dashed arrows illustrate trajectories of type I 
and type II. 



which translates into the vanishing of the curvature dm/ dx of the surface profile. The height profile represented by 
such a trajectory is a periodic array of step bunches which is everywhere smooth. In contrast, particle trajectories with 
positive total energy (type II solutions) reach the reflecting "hard wall" at u — at a finite speed. This implies that 
the surface curvature dm/dx — (du/dx)/2m diverges as m — ► 0, such that the surface profile develops singularities of 
the type 

h(x) - h{x ) ~ (x - x ) 3/2 (37) 

near the bunch edges at x = xq. This is the well-known Pokrovsky-Talapov law which describes how the rounded 
parts of an equilibrium crystal shape join a flat facel^, and is reflected in the scaling of the terrace sizes at the edges 
of the bunch (see Sect IV E(l . Because of the singularities at the bunch edges, type II height profiles have finite support 
and cannot be continued to describe a periodic array of bunches. 

We now show that the two types of trajectories imply different scaling properties for the step bunches. For a 
trajectory of type I to satisfy the boundary conditions u(0) = u(L) — for L — > oo, it is necessary to set the particle 
energy £ = 0. The right turning point is then located at the value u* = (3£>/2Jo) 2 at which V(u*) = 0, and the 
period of the trajectory is readily seen to be of the order of L ~ \J~BC j\ Jq\. Bunches of arbitrary (large) lateral size 
L can therefore be accommodated only by treating Jo a s a free integration constant, and to let Jo ~ l/L — > for 
L — > oo. This is the procedure adopted in Ref.34 for a slightly different case. In the present context it implies the 
scaling relation Z m i n ~ \j\fuF ~ | Jq| ~ and since the number of steps in the bunch is of the order N ~ L/l m [ n , 
we find that Z m i n ~ TV" 1 / 2 . For future reference we record the full expressions for L and / m i n , which read 

L=(216S)^N^l, U = ^=(^) 1/4 ^. (38) 

For type II trajectories with £ > 0, it is possible to make the period arbitrarily large while keeping Jo fixed, simply 
by increasing the energy. For large bunches we then have u max S> u* , and the decreasing part of the potential becomes 
irrelevant. The equation for the profile reduces to the form 

C ± ( m ^l) = -\J \m, (39) 
dx \ dx J 

which was first studied by Nozieres^, and leads to the scaling law^^Mi i min ~ 7V~ 2 / 3 (see Sect IV El for a detailed 
derivation). Anticipating the results of the numerical solution of the discrete step equations in Sect lVIll it turns out 



9 



that this scaling law is in agreement with the numerical data for / m ; n , while the prediction (|38|l for type I trajectories 
is not. We conclude, therefore, that the type II trajectories of l|36[) . with singular behavior at u = 0, are the relevant 
solutions for the description of step bunches. This immediately raises the question of how the current Jo, which then 
no longer can be treated as an integration constant, should be determined. The answer will be given in the next 
subsection. 

We add a final remark comparing the two types of solutions. On purely mathematical grounds, at first sight the 
type I solutions may seem to be preferable, because they avoid the singularities at the bunch edges and allow to 
describe a periodic array of many bunches. This is in fact not the case. For a periodic type I profile, the end of one 
bunch (the point where m = u = 0) defines the beginning of the next. As, according to (|38(l . the bunches steepen 
with increasing size, the mean slope of the surface also increases without bound. This is in contradiction to the time 
evolution of a real surface, for which the mean slope is fixed and the steepening of the bunches is compensated by the 
growth of large flat regions between the bunches (see Fig^). For type I solutions, the region between bunches shrinks 
to a point. Therefore they cannot be taken at face value as a global description of a surface with many bunches. Just 
like for the type II solutions, which terminate in singularities, also type I solutions have to be complemented by a 
separate description for the flat regions between bunches. 



D. The mean surface current 



Our strategy will be to fix Jo by analogy with two related problems, which are mathematically equivalent to the 
present one (see Sect I VII for further discussion). The two problems are the linear step growth model considered in 
Ref.42, and the model of surface electromigration in the attachment/detachment-limited regime considered in Ref.24. 
Indeed, the step equations derived by Liu and WeeksSi for surface electromigration in the presence of desorption but 
without Ehrlich-Schwoebel barriers are (apart from the different physical meaning of the coefficients) identical to our 
equations (TJJ. By appealing to the analogy with surface electromigration, we can associate a microscopic current ji 
with the terrace between the steps i and i + 1, which is given by the expression 



-h Q D[An(x i+ i) - A/i(xi)] 



(1-/3) 
2(1 + 0) 



h R e (x l 



(40) 



Here the abbreviation D = D s Q,n e s j3 /[ddk^T(\ + (3)] has been introduced. The physical meaning of this current is 
clear in the context of surface electromigration - it is simply the current of adatoms driven across the terrace by the 
combined action of the electromigration force and the chemical potential gradient. The notion of a surface current 
induced by attachment asymmetry is also well established in the context of epitaxial growthiM242i^24Si. It is less 
evident that the concept can be extended to sublimation, where the atoms detached from the steps do not necessarily 
remain on the terrace. However, as was noted above in Sect IV Al we are working here in a limit where the diffusion 
length X s is very large, and hence the mass transport is essentially confined to the surface. Mathematically, the 
derivation in Ref.42 shows that for a general set of step equations which arc linear in the terrace widths, 



at 



(41) 



the leading order (stabilizing or destabilizing) part of the surface current is proportional to the asymmetric combination 
7« — Id of the contributions from the two terraces, while the symmetric combination j u + -fd gives rise to the overall 
growth (or sublimation) rate of the surface, and to the symmetry-breaking part of the surface current. 
For a perfect step train with constant step spacing I, the current ijlU)) is equal to 



Jfla 



1-/3 
2(1 + /?) 



hoR e l, 



(42) 



which is just the first term in (|31[) evaluated at slope m = ho/l. Following Liu and Weeks^i, we now argue that 
the expression (|42|l remains valid also for a periodic array of step bunches. Each bunch contains ./V steps, and hence 
consists of N — 1 short terraces. The width of each bunch is L, and the bunches are separated by terraces of length 
L t . We denote by A/x_ and A/z+ the values of the chemical potential at the lower and upper edges of the bunch. 
Summing the expression (|40|l across the bunch, we find that the current in the bunch is equal to 



3b 



hp 
' N- 1 



(1-/3) 
2(1 + 13) 



R„L 



(43) 
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while the current on the terrace is 

j t = -haD{An- - A M+ ) - 2{l+% kokeLt - (44) 

Stationarity requires that jj, — jt, which yields an expression for the chemical potential difference A/i + — A/i„. 
Inserting this back into l|43|) or (|44|l . wc find the simple result 

• • (i-/ 3 ) u p L + L t (AKS 

Jb = Jt = -WW) hoRe ^^ J (45) 

which coincides with (|42|l . Thus the overall current remains constant, at its value for a regular step train, throughout 
the bunching process, and the appropriate expression to use for Jo in lj36|l or (|39|) is Jfl a t- 

It is possible to argue for this choice of Jo also without reference to the discrete step dynamics. Indeed, setting 
Jo = Jft&t ensures that the minimum of the potential V(u) is located at the value u = {ho /I) 2 corresponding to the 
mean surface slope. In this way the initial regular step train, which is clearly an (unstable) stationary solution of the 
discrete step dynamics, is retained as a solution also in the stationary continuum equation. 



E. Derivation of the scaling laws 

The solution of the mechanical problem (|36|l allows to express the bunch width L and the bunch height H in terms 
of the maximal slope m ma x = ^/it m ax- Using energy conservation and neglecting the linear term — Bu in the potential 
V(u), we obtain the expressions 




and 




The dimensionless coefficient in (|46|) was obtained by numerically evaluating the corresponding integral. Inserting 
the expressions Q29JI and (|42|l for C and Jo, respectively, we can write the minimal terrace length Z m j n and the bunch 
width L in terms of the number of steps and dimensionless ratios of length scales, as 

¥ = 21,5 (t^)" 3 (£)"* (t) 2 ' 3 (t) w ~ 2/s = m 



1 3.25 (-?-;) (-^) (^) f^)« 1/3 = 3.25 S l ' 3 N"\ u'l, 



and 



1-/3/ \ddj V 1 J \ l 

Equations (|48|l and (|49|) are the central results of this section, and will be compared to numerical simulations of the 
discrete step model in Sect. IVIll Together l|4*5|) and (|^9l imply the universal relation L/l nlin sa 1.29/V, independent 
of all physical parameters. This shows that the minimal terrace length is only a factor 0.78 smaller than the mean 
terrace length L/N within the bunch, hence most terraces in the bunch have a size of order l m in- 

A separate scaling law holds for the size of the first (and last) terrace at the edges of the bunch. To derive it, we need 
to determine the precise bunch profile near u = 0, which is of the general form (|37(l . Since the mechanical potential 
V(u) vanishes at u = 0, energy conservation implies that (C / A)(du/ dx) 2 = £ = V(u max ) for u — * 0. Integrating this 
from the bunch edge at x — yields the slope profile m(x) — (4£ /C) 1 / 4 ^ 1 / 2 , and correspondingly the height profile 
h(x) = (2/3)(4£/C) 1 / 4 a; 3 / 2 . There is an ambiguity in how to estimate the size l\ of the first terrace - a quantity 
manifestly related to the discreteness of the surface in the vertical direction - from these continuous profiles. Natural 
choices would be to require that (i) h{x = li) — ho, (ii) m{x = l\) = h^/li or (iii) m(h = ho) = ho/h- It is easy 
to check that all three choices imply a scaling relation of the form ~ S' 1 / 3 A r ~ 1 / 3 , but with different numerical 
coefficients. Because the continuum equation is derived 4 ^ from a set of discrete equations for the terrace sizes (the 
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inverse slopes) as a function of the layer height, we claim (and confirm numerically in Sect lVIlfl that choice (iii) is the 
appropriate one. This results in 

7 = 41/3 (r^)" (t)" (£)" (t) = 4 " 3 sl " N ^ < 50 > 

The same kind of analysis can be carried out for the type I solutions discussed in Sect IV Cl One finds that l\/l ~ 5 1 / 4 
independent of the bunch size N . 



F. General step interaction 

The above considerations are easily generalized to different values for the exponent n describing the decay of the 
step-step interaction as l/l n . The expression JHJ for the chemical potential at step i then becomes 



A//(xj) 
fc R T 



lo 



•^2 + 1 *^2 



n+l 



lo 



*^2— 1 



n+l 



(51) 



with l = {nQ.A/k^T) 1 / {n+1 \ and going through the manipulations of Sect II VI one obtains the generalized relaxation 
equation 



dh d_ 

dt dx 



with 



fdh\ 


-l 


" d (dh \ n ~ 


1 d 2 h 


\dx) 




dx \dx J 


dx 2 


(rH 


-iK 


+1 nD s n e s (3hl- 


-n 



dd(i + 0) 

The analysis of the resulting full continuum equation leads to the scaling relations^ 

for type II solutions. Specifically, the relations I|48I5U|) generalize to 

lmln/l = (165 n ) 1/( " +1) ^~ 2/(n+1) = (165„/iV 2 ) 1/(n+1) , h/l = (45„/iV) 1/( " +1) , 
where the dimensionless parameter S n is defined by 



S n — 



13 \il 



2;n+l 







K d K u t fl 



1 - (3 dl n + 2 K d -K u l \ l 



n+l 



(52) 



(53) 



(54) 



(55) 



(56) 



Comparing the two expressions in (|55fl it is seen that ^ min becomes smaller than Z x only for N > 4; this can be viewed 
as a necessary condition for the onset of scaling. If we require in addition that h/l < 1, it follows that such small 
bunches can form only provided S n < 1. 

For type I solutions the scaling relations corresponding to l|54|l read 



L ~ S 11 /^ 2 ) /V™/("+ 2 ) / 



S l/(n+2) N -2/(n+2)_ 



(57) 



G. Corrections to the asymptotic behavior 



The scaling laws derived in the preceding subsection are valid asymptotically for large bunches, when the two 
approximations made in their derivation are well justified: First, neglecting the symmetry-breaking contribution to 
the current (|3 1 1) . and second, neglecting the constant term in the particle equation Ij^tjl) . The second approximation is 
valid provided u max 3> u*, where u* = (3B/2 Jo) 2 is the value of u at which V(u*) = 0. Inserting the expression 142|) 
for Jo, we see that the corresponding slope m* = y/u* is simply of the order of the mean slope ho /I of the surface. 
Thus the linear term in V(u) can be ignored throughout the bunch. 
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FIG. 3: Stationary bunch shapes computed by numerical integration of Eg. 1581 with S = 0.14 and j3 = 0.01. The figure shows 
the dimensionless slope (l/ho)m(x) as a function of the rescaled coordinate x/l for different values of the maximum slope. The 
horizontal dotted line shows the value m = ho /I corresponding to the mean slope of the unperturbed surface. 



Including the symmetry-breaking term, the stationarity condition J — Jq can be brought into the form 

n d 2 v I. „ 1 + 1 dv , ras 
V^3 (1 -^ )+ 18(T-^^^ (58) 

where v = (l/ho) 2 u is a dimensionless version of u, normalized such that the mean slope corresponds to v = 1, and 
y = x/l is the dimensionless spatial coordinate. In the mechanical analog of Sect IV (^1 the symmetry-breaking term 
corresponds to a friction force which, because of its sign, acts "backward" in time. The friction force is proportional 
to v -3 / 2 . Hence it is negligible deep inside the bunch, where v ^> 1, but becomes important near the edges of the 
bunch. 

It is instructive to follow the solution of Ea. (|58|) starting from the center of the bunch, where v(y) takes its maximum 
value i> m ax and dv/dy = 0. Integrating forward in "time" y, towards the upper edge of the bunch, the friction force 
adds to the acceleration of the particle towards v = 0, which is therefore reached earlier than in the absence of the 
symmetry-breaking term. Moreover the singular behavior at v = is altered: Balancing the symmetry-breaking term 
against the intertial term on the left hand side of lf5~%|) , it is straightforward to show that the standard behavior i|57|) 
of the height profile is modified into 

h(x )-h(x) ~ (x -x)^ 3 . (59) 

Conversely, when moving backward in time (towards the lower edge of the bunch) the particle is delayed by the friction 
force. Since the friction coefficient diverges as v — > 0, the point v = is never reached. Instead, the trajectory bounces 
back and approaches the stable potential minimum at v — 1 in an overdamped or damped oscillatory manner. Since 
the parameter S plays the role of the particle mass in l|58|) , the effects of friction increase with decreasing S, while 
they decrease when increasing the initial particle energy (i.e., the value of v max , and, hence, the size of the bunch). 

Thus the symmetry-breaking term modifies the nature of the solutions of (|58l) in a qualitative way: Whereas the 
relevant solutions of the frictionless particle problem (|3GH have finite support in x (the trajectory returns to u = 
in a finite time), the solutions of l|58() extend all the way to y — — oo, where v attains the limiting value v = 1. In 
physical terms this implies that the bunch width L can no longer be sharply defined within the continuum theory. 
Nevertheless, the numerical solutions of (|58|) depicted in Fig. [31 show that this effect is completely negligible already 
for moderately large bunches and physically relevant values of S. Deviations from the solution of the frictionless 
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equation l|36|) occur only in the range v <C 1, which is irrelevant for the description of actual bunches (recall that 
v = 1 corresponds to the mean slope ho /I of the unperturbed surface). This conclusion is supported by a scaling 
analysis in the spirit of Ref.36 (see Sect lVIIl)) . 

We will see in Sect IVIll that the left-right symmetry of the bunches is indeed broken in a way that is qualitatively 
reminiscent of the solutions of (|58|l with very small S. However, the sign of the observed symmetry-breaking is 
opposite to that predicted by Q58JI. and we will argue that its origin is in fact completely different. 



VI. EQUIVALENT PROBLEMS 



We noted already in Sect lVDl that the equations of step motion (|10I15I16[1 are mathematically equivalent to ap- 
propriate limiting cases of those obtained for step bunching instabilities induced by electromigration and growth in 
the presence of inverse Ehrlich-Schwoebel barriers. Here we elaborate on that observation and translate the results 
derived from the continuum theory to the different physical contexts. 



A. Electromigration 

Discrete and continuum equations for electromigration-induced step bunching in the attachment /detachment limited 
regime have been derived by Liu and Weeks 2 ^, and their equations are readily seen to be of the same form as ours. In 
our setting, an electromigration force F acting on the adatoms can simply be added to the chemical potential gradient 
on the right hand side of (|24|) . This gives rise to an additional contribution Jp = aF to the surface current J in i|30[) . 
which is also inversely proportional to the slope m [see the expression <|25[) for a in the attachment/detachment limited 
case], and which is destabilizing for F < (force in the downhill direction). In the absence of an Ehrlich-Schwoebel 
effect (j3 = 1), the electromigration current is the only destabilizing contribution. The results of Sect Ivl carry directly 
over to this case, once the dimensionless parameter S has been identified along the lines of Sect IV Bl One finds the 
simple result 

(60) 

and hence from (|48I49I50() we obtain the predictions 

Un=f-^-J iV- 2 / 3 , L«3.25 i^—j N / , *!=(-—) N / . (61) 
Similar formulae have been reported previously in the literature. Sato and Uwaha derived the results 26 

1/3 /on A\ 1 / 3 /OO/lX 1 / 3 



L/N « 2.59 j N-^, U = j N-^ h = J N^, 



(62) 



which arc of the same form as the expressions in l|61|l . with the kinetic length d — d u — dd replacing the mean step 
spacing /. This reflects the fact that Sato and Uwaha work in the diffusion-limited regime d < I. 

Stoyanov and Tonchev 2 ^ have developed a continuum description for electromigration-induced step bunching in the 
diffusion-limited regime. Assuming the relation d = ay for the kinetic length, which holds for non-permeable steps in 
the absence of an additional barrier against attachment, they obtained the evolution equation 



dh d 
dt dx 



C d 



2 



= 0, (63) 



where m = dh/dx is the surface slope, and the coefficients B and C are given by 

- 2D a n e s Fna ± - 6D s n e s Q 2 A 

B = ' C= h i, r ■ ( 64 ) 

k B T h k B T 

The surface is unstable when the force is in the downhill direction, F < and B > 0. The relaxation term in 
(I63|l is simply the product of the chemical potential variation l|27|) and the expression l|26|l for the mobility in the 
diffusion- limited case, where it has been used that dd <C A s and /3 = 1. 
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The analysis of stationary solutions of (|63() proceeds along the lines of Sect IV Cl In fact, the mechanical analog 
resulting from setting the square bracket in (|63fl equal to a constant current Jo is formally identical to 1|36[) . with the 
potential V(u) = —Jqu + Bv?l 2 . Despite the formal similarity, however, the problem differs from that considered in 
Sect IV Cl in one important respect: Since the overall surface current is downhill, we have to chose Jo < 0, and hence 
both terms in the potential are positive; the potential has no minimum, and type I trajectories do not exist. For the 
type II trajectories the term proportional to Jo becomes irrelevant at large slopes, and hence one may as well set 
Jo = 0, as was done in Ref.29. The bunch shape is then given by the solution of l|39|) . and the results of Sect IV El can 
be taken over. This yields the scaling relations 

f A \ 3 fa A\ 3 

L/iV«2.58f^Lj TV" 2 / 3 , Z min = 2f-Lj at-2/3_ (65) 

Apart from a small difference in the dimensionless prefactor, the expression for L/N is identical to the one reported 
in Ref.29. 

Relations of the form (|65|l have been used in the interpretation of several experiments on silicon surfaces, where a 
scaling of the minimal terrace size^i and the mean terrace size 3 ^ as N~ 2 ^ 3 was observed. The similarity between the 
expressions <|61I62I65[) implies that it is not possible to distinguish between attachment/detachment limited kinetics 
and diffusion-limited kinetics on the basis of the observed scaling; see Sect I Villi for further elucidation of this point. 
However, the resulting estimate for the ratio A/F depends crucially on which kinetic regime is assumed. Consider 
for example the results obtained by Fujita, Ichikawa and StoyanoAS^i for the maximum bunch slope ho/l mm at 1250° 
C. Using the relation l|65() for the diffusion-limited regime yields the estimate F/A i=s 3 x lCU 6 nm- 2 ; this is somewhat 
larger than the value reported in Ref.31, because in that work the authors used an expression for the mean terrace 
width I = L/N. Because of the additional factor of I, application of (|61|) yields instead F/A «6x lO^nm -2 , which 
is smaller by a factor of 50. Correspondingly the effective valence Z* of the silicon adatoms, defined through the 
relation F = Z*eE between the electromigration force and the electric field E, will be smaller by the same factor. 



B. Growth with inverse Ehrlich-Schwoebel barriers 



Growth in the presence of an inverse Ehrlich-Schwoebel effect is described by the stationary diffusion equation © 
with R > and 1/t s — 0, and the boundary conditions © with K u > Kd, i.e f3 > 1. While the possibility of inverse 
ES barriers is debatable on the microscopic level, the inverse ES effect is may serve as a useful effective description 
of more complex step bunching mechanisms—. The equations of step motion can be found, e.g., in Ref.12. We 
consider the limiting case of fast attachment to the descending step and slow attachment to the ascending step, i.e. 
dd 3* Xi — Xi-i 3> d u . In this limit only the upper terrace contributes to the growth of the step, and the destabilizing 
part of the dynamics reduces to the linear form (|41|l with 7<j = and 7„ = —Rtt (note that in our setup the steps 
recede during growth and the upper terrace trails the step). The continuum equation is thus of the same form as in 
sublimation and electromigration, and the results of Sect fVl carry over with the identification 

of the dimensionless parameter. The application of the continuum theory to this problem will be the subject of a 
separate publication^. 



VII. NUMERICAL ANALYSIS OF THE DISCRETE STEP DYNAMICS 



Extensive numerical simulations of the discrete step dynamics have been carried out to test the predictions of the 
continuum theory. In this section we report on simulations for the sublimation problem; a comprehensive numerical 
study of growth in the presence of an inverse Ehrlich-Schwoebel effect will be presented elsewhere^. We work with 
cyclic boundary conditions, xm+i — x\ + Ml where I is the average intcrstcp distance of the vicinal surface and 
M is the total number of steps, and prepare the system in one of two kinds of initial conditions. Under natural 
bunching conditions, the integration is started from a vicinal surface with steps which slightly deviate from their 
regular positions. This leads to a surface consisting of many bunches of steps separated by large terraces, which 
slowly coarsens (FigQJ. On the other hand a single bunch can be prepared by chosing an initial step configuration, 
corresponding to a bunch, which contains almost all of the steps in the system. The integration then provides results 
for the steady state shape of the bunch and the average value of the number N of steps in it (the remaining M — N 
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FIG. 4: Profile of a crystal surface after some time of sublimation. The inset shows the enlargement of an individual bunch. 



steps are single - they are crossing the large terrace between the front edge of the bunch and its tail in our single 
bunch setup with cyclic boundary conditions). In both cases we assume that a given step belongs to the bunch when 
the distance to at least one of the neighbouring steps is smaller than 0.75 I. This definition is, in fact, arbitrary and 
it introduces some ambiguity in the results. The dependence of the minimum distance / m j n between the steps in the 
bunch on the number N of steps is less affected by this ambiguity than that of the total bunch width L, as in the 
former case only TV is influenced by the definition, whereas Z m ; n is precisely determined. 

We first present results from the numerical integration of the sublimation problem in case (b) of Sect lllll i.e. I|1(J|) 
was used with v u and Vd given by I|15I16|) , with the destabilizing terms depending linearly on the terrace widths. The 
equations of step motion contain 4 dimensionlcss parameters: X s /l, dd/l, lo/l and (3 = K u /Kd- It is reasonable to 
briefly discuss the values of these parameters. For instance the value X s /l = 100 means that A s = 1/zm at / = 10 ran. 
As far as the values of the parameter dd/l are concerned, it is difficult^ to evaluate the kinetic length dd- We should 
have in mind, however, that the Eqs. I|15I16|) are valid under the assumption [xi — iCi—i) <C dd and, therefore, we have 
to take dd/l ^> 1. The value of the parameter lo/l was assumed to be lo/l = 0.24 in order to keep the interstep distance 
in the bunch to be in a convenient interval. Finally, the values of (3 used were (3 = 0.01 and /3 = 0.1 corresponding 
to a rather high Ehrlich-Schwoebel barrier. Figure |S] shows four sets of data for l m i n as a function of N obtained 
for two different step interaction laws, n = 2 and n — 3, using the natural bunching as well as the single bunch 
initial conditions. In all cases excellent agreement with the theoretical prediction l|55l) is found. The same quality of 
agreement has been obtained for the problem of growth with inverse Ehrlich-Schwoebel barriers^. 

In Fig|H|we show data for the dependence of the total bunch width L on the number of steps. Although the overall 
magnitude of L is consistent with the prediction 149fl for type II stationary profiles (and rules out type I behavior), 
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FIG. 5: Numerical data for the minimum interstep distance, measured in units of I, as a function of bunch size. Open symbols 
show results obtained in the natural bunching geometry with 500 steps, while asterisks show data obtained from computations 
with a single bunch. The interaction strength was lo/l = 0.24 in all cases. Open squares and diamonds show data for dd/l = 10 
and X a /l = 100, triangles show data for dd/l = 150 and \ s /l = 200, and asterisks show data for dd/l = X s /l = 100; other 
parameters are given in the figure. Bold lines show the theoretical prediction Ij55^ for i m i n . 



a power law fit to the data yields an estimate a rs 0.44 which is intermediate between the type II (a = 1/3) and 
type I (a — 1/2) predictions. To gain some insight into this discrepancy, we take a closer look at the shapes of the 
bunches in the numerical simulation (Fig^J. It is clear that the bunches are distinctly asymmetric: While there is an 
abrupt change in the terrace length at the lower (trailing) edge of the bunch, at the upper (leading) edge the terrace 
lengths increase gradually. The asymmetry can be quantified by looking at the scaling of the size of the first (h) and 
last (Zjy) terrace in the bunch with the bunch size N (Fig[7|). While the data for l\ are in good agreement with the 
theoretical prediction l|55(l . the size of the last terrace is found to be essentially independent of N . Incidentally, the 
latter behavior is also characteristic of the type I stationary profiles (see Sect IV El) . More significantly, a constant last 
terrace size ijy ~ I results trivially from our way of numerically locating the bunch edge, if the terrace size increases 
continuously across the mean terrace size as one moves out of the bunch in the forward direction, i.e., if a sharply 
defined bunch edge in fact does not exist. 

What is the origin of the asymmetry in the bunch shape? We have shown in Sect IV (Jl that the symmetry-breaking 
term causes the bunch edge to fray out at one side, in a qualitatively similar manner to the behavior seen in Fig0| 
However, the blurring of the bunch edge is predicted to occur at the lower (trailing) edge, rather than at the upper 
edge, and in addition the effect becomes negligibly small already for moderately sized bunches (see FigEJ). We believe, 
instead, that the bunch asymmetry is intimately related to the exchange of crossing steps between bunches. These 
steps gradually accelerate out of the bunch at the leading edge, which translates into a gradual increase of the mean 
terrace size. Conversely, when approaching the next bunch from behind, the crossing step decelerates quite abruptly, 
because it is primarily fed from behind (this is particularly true in the almost one-sided regime mainly considered 
in our simulations) . The exchange of steps between bunches also implies that the bunches move laterally at a speed 
which is different from the mean sublimation rate. To capture the asymmetry in the continuum theory, it will therefore 
be necessary to go beyond the stationary solutions considered in Sect0 which by construction are symmetric, and 
to investigate solutions describing moving and interacting bunches. It is worth pointing out that the existence of 
crossing steps partly invalidates the argument used in Sect IV E)l to fix the mean surface current, because the argument 
assumes that all steps reside in bunches^*. 

We close this section with some remarks concerning the limiting case (a) of Section ITTT1 We have shown that no 
bunching occurs in this case, i. e., when the assumption Xi — Xi-\ > dd is fulfilled, if only the linear and quadratic 
terms [(xi — Xi-\)/\ s and (x^ — Xi-i) 2 /\ 2 S ] in the expressions for the step velocity are taken into account. To clarify 
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FIG. 6: Numerical results for the bunch width L as a function of bunch size. Data are shown for two of the parameter sets 
displayed in FigEl P = 0.1, \ s /l = 200, d/l = 150 (diamonds), and f3 = 0.01, \ s /l = 100, d/l = 10 (crosses). In both cases 
lo/l = 0.24 and the system contained 500 steps. Dashed and full lines show the predictions for type I and type II solutions, 
respectively. A power law fit to the data yields L ~ iV ' 44 . 



the problem of step bunching instability in the limiting case (a) we did a lot of numerical work making use of the full 
expressions for the step velocity. Integration of the equations of step motion proved the existence of step bunching 
at parameter values (3 = 0.01, lo/l — 0.003 and dd/l = 1/3, dd/l — 1/30 and dd/l = 1/300. It is essential to note, 
however, that the magnitude of the step-step repulsion energy used in these integration runs was much smaller than 
in the integration of the equations obtained in the case dd ^> I, where we used lo/l = 0.24. On the other hand, the 
density of steps in the bunch is higher in the case dd 3> I compared with the opposite limiting case dd <SC I- These 
findings indicate a strong impact of the parameter dd on the bunching process. When the parameter dd is larger than 
the interstep distance I, bunching occurs even at a very strong repulsion between the steps and the step density in 
the bunch is rather high. On the contrary, when the parameter dd is smaller than the interstep distance I, bunching 
occurs only at a very weak repulsion between the steps and the step density in the bunch is relatively small. This is 
not surprising because in the case under consideration neither the linear nor the quadratic term induces instability 
of the vicinal surface. Destablising terms are of higher order, i. e. [(xi — £i_i)/A s ] 1 ' with v > 2 and their effect is 
relatively weak so that it cannot dominate a strong repulsion between the steps. It is interesting to note, however, 
that in this case of weak instability (dd <C I) the minimum interstep distance in the bunch scales with the number of 
steps in exactly the same way (l m i n ~ N~ 2 / 3 ) as in the case of strong instability (dd 3> I). Further discussion of this 
regime will be presented elsewhere. 



VIII. UNIVERSALITY CLASSES OF STEP BUNCHING 



Before drawing some general conclusions in the next section, here we wish to put our work into the context of 
a classification scheme for step bunching instabilities proposed by Pimpinclli, Tonchev, Videcoq and Vladimirova 
(PTVV)i2ii. It is based on a generic continuum equation of the form 



dh 



d_ 

dx 



Kim s + ^§^m r 
mr ox z 



— const. 



(67) 



Here K\ and K 2 are material constants with K 2 > and Kiq > 0, the slope m = dh/dx is assumed positive 
everywhere, and g, k and n are exponents characterizing a class of step bunching instabilities. The exponent n is 
simply the exponent of the repulsive step interaction, and the exponent k reflects the slope dependence of the surface 
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FIG. 7: Scaling of the first and last terrace size in the bunch with bunch size for the same parameter sets as in Fig|S| The full 
lines show the theoretical prediction for type II solutions. 



mobility; k = 1 and k = correspond to slow and fast detachment/attachment kinetics, respectively (see Sect llVII . 
Equation (j67(l is a slight generalization 39 of the equation proposed in Ref.36, where only the case k = Q was considered. 

PTVV argued that the characteristic scaling exponents a and z introduced in Sect[I]can be extracted from l|67l) by 
requiring that the equation should be invariant under the scale transformation 



h(x,t) => b- a h{bx,b z t) 



for an arbitrary scale factor b. This yields the expressions 



7 = 



2(1 + n - k - 2p) 



(68) 



(69) 



where the scaling relation © has been used. 

Apart from the symmetry- breaking term in l|31(l . the continuum equation (j30|) for the sublimation problem is of 
the generic form l|67|l with g = — 1 and k = 1. It is straightforward to check that, under the rescaling (|68|l . the 
symmetry-breaking term is smaller by a factor of b~ a compared to the leading term ~ 1/m, and hence it becomes 
negligible at large scales; this is consistent with the detailed analysis in Sect IV Cl For g = — 1 and k = 1 the exponents 
in 169|) reduce to a = (n + 2) /n and 7 = 2/(2 + n), which we readily recognize as the scaling exponents characteristic 
of type I solutions [compare to i|57|) ]. To obtain the exponents for type II solutions (which, as was shown in Sect lVIll 
correctly describe the bunch shape), we have to set g — instead of g = — 1 in i|rj9"|) . 

The reason for the shift in the value of g is evident in view of the considerations of Sect IV Cl The scaling argument 
of PTVV assumes that all terms in (|67|l are of a similar order of magnitude; in particular, the stabilizing term is 
balanced against the destabilizing current Bm e . However, for type II solutions the total current is fixed at a value Jo 
which is independent of the slope, and which dominates over the term Bm e for large slopes when g < 0. Thus the 
stabilizing term is balanced against a slope-independent current, which effectively implies that g = 0. 

The argument clearly extends to any negative value of g, and suggests that generally g should be replaced by 
max(0, g). For the continuum evolution equation i|63|) describing electromigration- induced step bunching in the 
diffusion-limited regime, which corresponds to g = 1 and k = 0, the ambiguity regarding type I and type II solutions 

Since the static exponents a and 7 in (|69|l only depend on the sum g - 



does not arise (see Sect lVI All . Since the static exponents a and 7 in (|69|l only depend on the sum g + k, it is evident 
that they take the same values for g = 0, k = 1 as for g = 1, k = 0; for this reason the static scaling exponents for 
electromigration- induced step bunching take the same values in the diffusion- limited and the attachment/detachmcnt- 
limited regimes. 

The scaling theory of PTVV also makes predictions about the coarsening behavior of the bunched surface, which 
our analysis of stationary bunch shapes clearly cannot address. Setting g = — 1, the expressions (|69|l yield a/z = 1/2 
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FIG. 8: Time dependence of the mean bunch size for the parameter sets in Figs|S|and|7| and an additional set with /3 = 0.1, 
As/l = 100, d/l = 100, Iq/1 — 0.24, and step interaction exponent n = 3. The bold lines illustrate the predictions of the scaling 
theory for g — — 1 {a/z = 1/2) and for g = 0, n = 2 (a/z = 3/4). 

for the exponent in the coarsening law Q, which is independent of both n and fe; on the other hand, with g = one 
obtains a/z = (n + 1) /2n for fc = 1. In FiglHlwe compare numerical data for the temporal evolution of the mean bunch 
size to these two coarsening laws. The simulations seem consistent with the "superuniversal" value a/z — 1/2, but 
also a/z = 3/4 (for n — 2) or 2/3 (for n = 3) cannot be ruled out. More extensive simulations arc needed to firmly 
pin down the coarsening behavior; this is particularly true here because, in contrast to the static scaling properties 
discussed earlier in Sect lVIll we do not have any analytic information about the coefficient of the coarsening law. In 
a recent study of a simple toy model of step bunching, which ignores the repulsive step-step interactions and allows 
steps to coalesce, it was necessary to go to extremely long times, equivalent to the growth of more than 10 5 monolayer, 
to ascertain the true asymptotic coarsening behavioi^,. 



In this paper we have presented a detailed analysis of the step bunching instability caused by an Ehrlich-Schwoebel 
effect during sublimation in the limit of a small desorption rate. When the kinetic length dd = D s /Kd is large 
compared to the average distance between the steps, the instability is strong and bunches of steps appear even at 
strong repulsion between the steps. In the opposite case (kinetic length dd smaller than the interstep distance) the 
instability is weak and bunches occur only when the step-step repulsion is several orders of magnitude weaker than 
in the previous case. 

A central part of the work is the derivation of the continuum evolution equation in Section^ and the careful analysis 
of its stationary bunch solutions. We have shown that two different types of stationary solutions with different scaling 
properties can be found, depending on whether the mean surface current Jo is kept fixed or not. Following Ref.24, we 
have argued that J is independent of bunch size, and that the correct bunch shape is given by the type II solution, 
which describes a bunch of finite extent with Pokrovsky-Talapov-type singularities at the edges. This is confirmed by 
the excellent agreement with numerical simulation results for the minimal interstep spacing Z m i„ and the first interstep 
spacing in the bunch lx presented in Sect lVIll 

On the other hand, we find noticeable deviations of the behavior of the total bunch width L from the type II 
prediction. We suggest that the discrepancy may be related to the distinct asymmetry between the leading and the 
trailing edges of the bunch: The terraces between the crossing steps escaping from the leading edge of the bunch 
appear to contribute strongly to the total bunch width, to the extent that asymptotically L may be considerably 
larger than Nl m i n . Further clarification of the issue requires a better understanding of the motion of bunches and the 
interactions between bunches, which is beyond the scope of the present paper. 

Our work has important consequences for the recently proposed scenario of universality classes for step bunching 
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instabilities 36 . First, we have pointed out the mathematical equivalence between appropriate limits of the step 
bunching instabilities caused by sublimation, growth and electromigration on the levels of both discrete step dynamics 
and continuum evolution equations. This equivalence gives a very clear meaning to the notion of a universality class, 
and we believe that the particular class considered in this paper (characterized, in essence, by the linearity of the 
destabilizing terms in the step equations) is in fact relevant to a wide range of experimentally realized systems. As 
a practical matter, the unified view provided by the continuum approach allows us to derive explicit formulae for 
the bunch shape which, through the identification of the scaling parameter S, are directly applicable to these diverse 
realizations of step bunching. Second, we have shown that the procedure employed in Ref.36 to extract the scaling 
exponents from the continuum evolution equation captures only one type of solution (the type I solutions of Sect IV Cjl . 
which is not the relevant one at least as far as the time-independent scaling properties are concerned. 

A crucial question that should be addressed in future work concerns the coarsening behavior of the bunched sur- 
face, and the relationship between coarsening dynamics and bunch motion. As was discussed in Sect lVIlil the 
present work remains inconclusive on this point. It is remarkable, however, that a very robust scaling of the 
mean bunch size and bunch spacing as N ~ £ ~ i 1 / 2 has been observed in a number of numerical simulations, 
both for electromigrationSSiSiSii and growth with inverse ES barriersi 2 *^, as well as in an experimental study of 
electromigration-induced step bunching on Si(lll)iS. Liu and Weeks^ have proposed an elegant explanation for the 
ubiquity of the t^-scaling within a continuum setting; their argument presupposes, however, as does the scaling 
approach of PTVV^ 6 ., that the bunch spacing £ is the only lateral length scale in the problem, although the bunch 
width L clearly comes into play as wells*. This remains true even if the internal bunch structure is eliminated by 
allowing the steps to coalesce^. Thus the origin of the observed temporal scaling remains to be understood. 
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